function [XBoot, RmBoot, RfBoot] = BootGenData(X, Rm, Rf, NBoot)
% DGP is assumed to be iid for X, Rm, and Rf

[T, K] = size(X);
RmBoot = NaN(T, NBoot); RfBoot = NaN(T, NBoot); XBoot = NaN(T, NBoot, K);

% If X is a matrix, then it can happen that it is not a full panel.
% Meaning some values X(t,k) are missing. The code handles this.
% And generates the bootstrapped sample that is also NOT full panel.

%% Generate in-sample estimates
muRm = mean(Rm); eRm = Rm - muRm;
muRf = mean(Rf); eRf = Rf - muRf;
muX = mean(X,'omitnan'); eX = X - muX;

%% Generate bootstrapped data
Textra = ceil(T/10); % these are the 10% extra observations to generate
% Now generate NBoot*(T+Textra) random numbers between 2 and T
TIndex = randi([2 T],NBoot*(T+Textra),1); TIndex = reshape(TIndex,T+Textra,NBoot);
for b = 1:NBoot
    Tb = TIndex(:,b);
    eRmb = eRm(Tb); eRfb = eRf(Tb); eXb = eX(Tb,:);
    Rmb = muRm + eRmb; RmBoot(:,b) = Rmb(Textra+1:end);
    Rfb = muRf + eRfb; RfBoot(:,b) = Rfb(Textra+1:end);
    for k = 1:K
        Xb = muX(k) + eXb(:,k); XBoot(:,b,k) = Xb(Textra+1:end);
    end
end
return